#!/usr/bin/env python
"""

    parse_hmm_search_result.py

"""

################################################################################
#
#   parse_hmm_search_result
#
#
#   Copyright (C) 2009 Leo Goodstadt
#
#   This program is free software; you can redistribute it and/or
#   modify it under the terms of the GNU General Public License
#   as published by the Free Software Foundation; version 2
#   of the License
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#################################################################################

import sys, os
# add self to search path for testing
if __name__ == '__main__':
    exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
    sys.path.append(os.path.abspath(os.path.join(exe_path,"..", "python_modules")))
    module_name = os.path.split(sys.argv[0])[1]
    module_name = os.path.splitext(module_name)[0];
else:
    module_name = __name__

    
from collections import defaultdict
import re
import json
from nested_dict import nested_dict

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions        


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# setup globals
domain_scores_header_re = []
domain_scores_header_re.append(re.compile(r"Sequence +Domain +seq-f +seq-t +hmm\-f +hmm\-t +score +E\-value"))
domain_scores_header_re.append(re.compile(r"(?:\-+ +){7}\-+"))
seq_name_re = re.compile(r"([^:]+):(\d+):(\d+):frame=([012]):strand=([01])")
query_hmm_str = "Query HMM:"
query_hmm_str_len  = len(query_hmm_str)


#_________________________________________________________________________________________    

#   _check_alignment_id

#_________________________________________________________________________________________    
def _check_alignment_id (hmm_results, current_hmm, line, cnt_line):
    # handle Id line
    id_fields = line.split()
    if id_fields[5] != "from" or id_fields[7] != "to":
        raise Exception("Unexpected file format on #%i\n[%s]\n" 
                                % (cnt_line - 3, line.rstrip("\n")))
    beg, end = int(id_fields[6]), int(id_fields[8][:-1])

    seq_name = id_fields[0][0:-1]

    if len(hmm_results[current_hmm][seq_name][beg, end]) > 7:
        raise Exception("hmm alignments already parsed for s %i->%i on #%i\n[%s]\n" 
                                % (seq_name, beg, end, cnt_line, line))
    
    if len(hmm_results[current_hmm][seq_name][beg, end]) != 7:
        raise Exception("No domain hits recorded for hmm_alignment matching %s:%i-%i on #%i\n[%s]\n" 
                                % (seq_name, beg, end, cnt_line, line))
    return seq_name, (beg, end)

#_________________________________________________________________________________________    

#   _handle_current_alignment

#_________________________________________________________________________________________    
def _handle_current_alignment(hmm_results, current_hmm, curr_seq_name, curr_beg_end, saved_hmm_match, saved_seq_match, cnt_line):
    """
    Save alignments
    """
    # make sure somewhere to save to
    if curr_seq_name == None:
        if len(saved_hmm_match) or len(saved_seq_match):
            raise Exception("Unexpected file format\n  "             +
                            "hmm_alignment = %s\n" % saved_hmm_match +
                            "seq_alignment = %s"   % saved_seq_match)
        #
        #  first alignment: nothing to save
        return
    
    # handle HMM line
    if saved_hmm_match[0:3] != "*->" or saved_hmm_match[-3:] != "<-*":
        raise Exception("Unexpected format for hmm alignment %s:%i-%i on line#%d [%s]\n" 
                                % (curr_seq_name, curr_beg_end[0], curr_beg_end[1],
                                    cnt_line, saved_hmm_match))
    saved_hmm_match = saved_hmm_match[3:-3]

    if len(saved_hmm_match) != len(saved_hmm_match):
        raise Exception("Alignment lengths do not match for %s:%i-%i\n%s ~ \n%s\n" 
                                % (curr_seq_name, curr_beg_end[0], curr_beg_end[1],
                                    saved_hmm_match, saved_seq_match))
    hmm_results[current_hmm][curr_seq_name][curr_beg_end].extend((saved_hmm_match, saved_seq_match))


#_________________________________________________________________________________________    

#   _handle_domain_alignments

#_________________________________________________________________________________________    
def _handle_domain_alignments(stream, hmm_results, current_hmm, cnt_line):
    """
    Handle part of HMM results which look like this:
        Alignments of top-scoring domains:
        chr10_random:0:113275:frame=1:strand=1: domain 1 of 1, from 11388 to 11859: score -379.3, E = 6
                           *->kaearsakntvdsYapylvsCpsddinlvRsamalsdnEseflqk..
                              +  a+ ++ +++s     v+C   +      a+     +++ ++k  
          chr10_rand 11388    LXNASPTFLPQPSSLXVRVTC---Y---WXIATVETKDQKKHRKKfy 11428
        
                           faTvdnntdddcptfigCvacAiirRklealkknvtlpsecsqcFadYc<
                              + ++    c   +  ++ A+       l  + ++p  c+  + +   
          chr10_rand 11814 LQKRVSSQRXGC--DTYLILTAVWLCQACCL-EGRSRPGHCEHIIXQIA  11859
        
                           -*
                             
          chr10_rand     -    -    

        chr10_random:0:113275:frame=1:strand=1: domain 1 of 1, from 11388 to 11859: score -379.3, E = 6
                           *->kaearsakntvdsYapylvsCpsddinlvRsamalsdnEseflqk..
                              +  a+ ++ +++s     v+C   +      a+     +++ ++k  
          chr10_rand 11388    LXNASPTFLPQPSSLXVRVTC---Y---WXIATVETKDQKKHRKKfy 11428
        
                           figCvacAiirRklealkknvtlpsecsqcFadYc<-*
                           +++        Rk e+      ++++    F+ Y    
          chr10_rand  4019 WLVTPQESHKLRKPENG--SCETSEKFKLLFQSYN    4051 
          

  Histogram of all scores:
    """
    # 
    # handle alignments to HMM hits
    # 
    LINE_TYPE_HMM, LINE_TYPE_MATCH, LINE_TYPE_SEQ = range(3)
    current_line_type = LINE_TYPE_HMM
    saved_hmm_match   = ""
    saved_seq_match   = ""
    curr_seq_name     = None
    curr_beg_end      = None
    curr_match_id     = None
    curr_hmm_id       = None
    for line in stream:
        cnt_line += 1

        line = line.rstrip("\n")
        
        # ignore empty line
        if not len(line):
            continue

        #
        #   end of alignment
        #   
        if line[0] != " ":
            _handle_current_alignment(hmm_results, current_hmm, curr_seq_name, curr_beg_end, saved_hmm_match, saved_seq_match, cnt_line)

            
            # 
            # end of whole section            
            # 
            if line == "Histogram of all scores:":
                return cnt_line
                
            # 
            # end of whole section            
            # 
            if "[no hits above thresholds]" in line:
                return cnt_line
            #
            # Make sure that the sequence identifier / positions match those previously
            #   parsed in the domain hit list
            # 
            if current_line_type != LINE_TYPE_HMM:
                print >>sys.stderr, "!!!!!!\n", current_line_type, cnt_line, line
            current_line_type = LINE_TYPE_HMM
            curr_seq_name, curr_beg_end = _check_alignment_id (hmm_results, current_hmm, line, cnt_line)
            saved_hmm_match = ""
            saved_seq_match = ""
            continue
        

        line = line.strip()

        if current_line_type == LINE_TYPE_HMM:
            saved_hmm_match += line
            
        elif current_line_type == LINE_TYPE_SEQ:
            fields = line.split()

            # save matching sequence
            if len(fields) == 4:
                abbrev, ignore_beg, seq, ignore_end = fields
                abbrev_len = len(abbrev)
                if abbrev != curr_seq_name[0:abbrev_len]:
                    raise Exception("Abbreviated sequence name does not match full name " +
                                    "on #%i\n[%s]\n" % (cnt_line, line))
                    
                saved_seq_match += seq

            # Last line of matching sequence might be empty to accomodate special
            #   characters on hmm match line
            elif len(fields) == 3:
                if fields [1] + fields[2] != "--":
                    raise Exception("Unexpected file format on #%i\n[%s]\n" % (cnt_line, line))
            else:
                raise Exception("Unexpected file format on #%i\n[%s]\n" % (cnt_line, line))
            

        # prepare for next line type
        current_line_type = (current_line_type + 1) % 3
            

    return cnt_line
    

#_________________________________________________________________________________________    

#   _handle_domain_scores

#_________________________________________________________________________________________    
def _handle_domain_scores (stream, hmm_results, current_hmm, cnt_line):
    """
    Handle part of HMM results which look like this:
        Parsed for domains:
        Sequence Domain  seq-f seq-t    hmm-f hmm-t      score  E-value
        -------- ------- ----- -----    ----- -----      -----  -------
        prot1    3/17    13380 13402 ..     1    23 []    28.7  1.4e-08
        prot2    5/19    13380 13402 ..     1    23 []    28.7  1.4e-08
        [no hits above thresholds]    
    """
    # 
    # make sure results header match
    # 
    for i in 0,1:
        line = stream.next()
        cnt_line += 1
        if not domain_scores_header_re[i].search(line):
            raise Exception("Unexpected file format on #%i\n[%s]\n" 
                                    % (cnt_line, line.rstrip("\n")))

    # 
    # handle HMM hits
    # 
    for line in stream:
        cnt_line += 1

        # hits end on an empty line
        if line == "\n":
            break

        #chr6_random:0:1875562:frame=0:strand=0   1/1   171276 171300 ..     1    23 []     6.4     0.07
        data = line.split()
        if len(data) != 10:
            if "no hits above thresholds" in line:
                break
                
            else:
                raise Exception("Unexpected file format on #%i\n[%s]\n" 
                                    % (cnt_line, line.rstrip("\n")))
        (   seq_name, domain_index, 
            beg, end, dots, 
            hmm_beg, hmm_end, 
            brackets, bit_score, e_value) = data
        beg, end = map(int,  (beg,  end))
            
        hmm_results[current_hmm][seq_name][beg, end]=[beg, end, hmm_beg, hmm_end, cnt_line, bit_score, e_value]
    return cnt_line
    
        
#_________________________________________________________________________________________    

#   parse_hmm_search_results

#_________________________________________________________________________________________    
from collections import namedtuple

t_hmm_result = namedtuple("t_hmm_result", "current_hmm, seq_name, beg, end, hmm_beg, hmm_end, cnt_line, score, e_value, aligned_hmm, aligned_seq")

def parse_hmm_search_results (stream):
    
    """
    Parse Results from running HMMer
    """
    stream = iter(stream)
    
    
    hmm_results = nested_dict(3,list)
    
    prev_pos = 0    
    current_hmm = "Unknown Query HMM"
    cnt_line = 0
    for line in stream:
        cnt_line += 1

        # current HMM. E.g.:
        #Query HMM:   smart_znf_c2h2_seed_alignment
        if line[0:query_hmm_str_len] == query_hmm_str:
            current_hmm = line[query_hmm_str_len:].strip()
            continue
            
        if line == "Parsed for domains:\n":
            cnt_line = _handle_domain_scores (stream, hmm_results, current_hmm, cnt_line)
        elif line == "Alignments of top-scoring domains:\n":
            cnt_line = _handle_domain_alignments(stream, hmm_results, current_hmm, cnt_line)
            for (current_hmm, seq_name, (beg, end)), (beg, end, 
                hmm_beg, hmm_end, cnt_line_, score, e_value, 
                aligned_hmm, aligned_seq) in hmm_results.iteritems_flat():
                yield t_hmm_result(current_hmm, seq_name, beg, end, hmm_beg, hmm_end, 
                                    cnt_line_, score, e_value, aligned_hmm, aligned_seq)
            hmm_results = nested_dict(3,list)
            




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Testing


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
import unittest
class Test_parse_hmm_search_result(unittest.TestCase):
    t_hmm_result_per_seq = namedtuple("t_hmm_result_per_seq", "beg, end, frame, cnt_line, score, e_value, aligned_hmm, aligned_seq")
    #_________________________________________________________________________________________    

    #   parse_seq_name

    #_________________________________________________________________________________________    
    seq_name_re = re.compile(r"([^:]+):(\d+):(\d+):frame=([012]):strand=([01])")
    @staticmethod
    def parse_seq_name(seq_name):
        """
        Get contig name, loci, frame and strand encoded in the sequence name
        """
        m = seq_name_re.match(seq_name)
        if not m:
            raise Exception("Unexpected sequence name format %s"  % seq_name)

        contig, contig_beg, contig_end, frame, strand = m.groups()

        return contig, contig_beg, contig_end, frame, strand

    @staticmethod
    def adjust_results_with_contig_offsets ( hmm_results):
        """
        If the translated sequence is broken up into fragments, calculate true offsets
        """
        adjusted_hmm_results = defaultdict(list)

        #debug_line_cnt = 0
        for (hmm_model, seq_name, beg, end, hmm_beg, hmm_end, 
            cnt_line, score, e_value, align_hmm, align_seq) in hmm_results:
            contig, contig_beg, contig_len, frame, strand = Test_parse_hmm_search_result.parse_seq_name(seq_name)
            contig_beg, contig_len, beg, end, hmm_beg, hmm_end, frame, strand = \
                    map(int, [contig_beg, contig_len, beg, end, hmm_beg, hmm_end, frame, strand])
            contig_end = contig_beg + contig_len
            beg -= 1
            hmm_beg -= 1
            if strand:
                m_beg = contig_beg + beg * 3 + frame
                m_end = contig_beg + end * 3 + frame
            else:
                m_beg = contig_end - end * 3 - frame
                m_end = contig_end - beg * 3 - frame
                print "mbeg = %d, contig_end = %d, end = %d "% (m_beg, contig_end, end)
            adjusted_hmm_results[hmm_model, contig, strand].append(
                                       Test_parse_hmm_search_result.t_hmm_result_per_seq(m_beg, m_end, frame, 
                                                    cnt_line, score, e_value, align_hmm, align_seq))


            #debug_line_cnt += 1
            #if debug_line_cnt > 10000:
            #    break

        #
        #   sort hits in strand order
        # 
        for (hmm_model, contig, strand) in adjusted_hmm_results.keys():
            adjusted_hmm_results[hmm_model, contig, strand].sort()
            
        return adjusted_hmm_results.items()


    #       self.assertEqual(self.seq, range(10))
    #       self.assert_(element in self.seq)
    #       self.assertRaises(ValueError, random.sample, self.seq, 20)


    
    def test_function(self):
        """
            test 
        """
        test_str = """hmmsearch - search a sequence database with a profile HMM
        HMMER 2.3.2 (Oct 2003)
        Copyright (C) 1992-2003 HHMI/Washington University School of Medicine
        Freely distributed under the GNU General Public License (GPL)
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        HMM file:                   /net/cpp-group/Leo/inprogress/6frame_translation_hmm/smart_znf_c2h2_seed_alignment.hmm [smart_znf_c2h2_seed_alignment]
        Sequence database:          /net/cpp-group/Leo/inprogress/6frame_translation_hmm/working_dir/hs_ncbi36_softmasked/hs_ncbi36_softmasked.00017.chr19.translated_fa
        per-sequence score cutoff:  [none]
        per-domain score cutoff:    [none]
        per-sequence Eval cutoff:   <= 10        
        per-domain Eval cutoff:     [none]
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        Query HMM:   smart_znf_c2h2_seed_alignment
        Accession:   [none]
        Description: [none]
          [No calibration for HMM; E-values are upper bounds]

        Scores for complete sequences (score includes all domains):
        Sequence                          Description           Score    E-value  N 
        --------                          -----------           -----    ------- ---
        chr19:0:63811651:frame=0:strand=0                     20471.0          0 884
        chr19:0:63811651:frame=1:strand=0                     20355.2          0 897
        chr19:0:63811651:frame=2:strand=1                     19683.6          0 856
        chr19:0:63811651:frame=0:strand=1                     19269.0          0 861
        chr19:0:63811651:frame=1:strand=1                     16303.5          0 759
        chr19:0:63811651:frame=2:strand=0                     13238.4          0 631

        Parsed for domains:
        Sequence                          Domain  seq-f seq-t    hmm-f hmm-t      score  E-value
        --------                          ------- ----- -----    ----- -----      -----  -------
        chr19:0:63811651:frame=1:strand=1 251/759 7857499 7857542 ..     1    23 []     2.0      1.2
        chr19:0:63811651:frame=1:strand=0 723/897 15040949 15040972 ..     1    23 []     2.0      1.2

        Alignments of top-scoring domains:
        chr19:0:63811651:frame=1:strand=1: domain 251 of 759, from 7857499 to 7857542: score 2.0, E = 1.2
                           *->ykCpigCg.....................ksfssk.aLkrHmrvH<-
                              ykC++ Cg+  +++++++ ++ +++++ +  f + ++L +H r H  
          chr19:0:63 7857499    YKCEE-CGkafnkspdhtxeklhkpercdNAFDNTsNLSKHKRNH   7857542

                           *
                           
          chr19:0:63     -   -    

        chr19:0:63811651:frame=1:strand=0: domain 723 of 897, from 15040949 to 15040972: score 2.0, E = 1.2
                           *->ykCpigCgksfssk.aLkrHmrvH<-*
                              + C + C+ +++ +++L  H  +H   
          chr19:0:63 15040949    LSCLHSCPACTPAPtHLHTHEHIH    15040972

        Histogram of all scores:
        score    obs    exp  (one = represents 1 sequences)
        -----    ---    ---
         -4      1      -|=                                                          
         -3      1      -|=                                                          
         -2      0      -|                                                           
         -1      2      -|==                                                         
          0      1      -|=                                                          
          1      0      -|                                                           
          2      1      -|=                                                          
        
        
        % No statistical fit available
        
        Total sequences searched: 6
        
        Whole sequence top hits:
        tophits_s report:
           Total hits:           6
           Satisfying E cutoff:  6
           Total memory:         21K
        
        Domain top hits:
        tophits_s report:
           Total hits:           6
           Satisfying E cutoff:  6
           Total memory:         22K
        """
        test_str = [s + "\n" for s in test_str.replace("\n        ", "\n").split("\n")]
        res = list(parse_hmm_search_results (test_str))
        self.assertEqual(res, [
                            t_hmm_result(   current_hmm='smart_znf_c2h2_seed_alignment', 
                                            seq_name='chr19:0:63811651:frame=1:strand=1', 
                                            beg=7857499, 
                                            end=7857542, 
                                            hmm_beg='1', 
                                            hmm_end='23', 
                                            cnt_line=32, 
                                            score='2.0', 
                                            e_value='1.2', 
                                            aligned_hmm='ykCpigCg.....................ksfssk.aLkrHmrvH', 
                                            aligned_seq='YKCEE-CGkafnkspdhtxeklhkpercdNAFDNTsNLSKHKRNH'), 
                            t_hmm_result(   current_hmm='smart_znf_c2h2_seed_alignment', 
                                            seq_name='chr19:0:63811651:frame=1:strand=0', 
                                            beg=15040949, 
                                            end=15040972, 
                                            hmm_beg='1', 
                                            hmm_end='23', 
                                            cnt_line=33, 
                                            score='2.0', 
                                            e_value='1.2', 
                                            aligned_hmm='ykCpigCgksfssk.aLkrHmrvH', 
                                            aligned_seq='LSCLHSCPACTPAPtHLHTHEHIH')])

    def test_function2(self):
        """
            test 
        """
        test_str = """hmmsearch - search a sequence database with a profile HMM
        HMMER 2.3.2 (Oct 2003)
        Copyright (C) 1992-2003 HHMI/Washington University School of Medicine
        Freely distributed under the GNU General Public License (GPL)
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        HMM file:                   /net/cpp-group/Leo/inprogress/6frame_translation_hmm/smart_znf_c2h2_seed_alignment.hmm [smart_znf_c2h2_seed_alignment]
        Sequence database:          /net/cpp-group/Leo/inprogress/6frame_translation_hmm/working_dir/mm_mm9_softmasked/mm_mm9_softmasked.00021.chr14.translated_fa
        per-sequence score cutoff:  [none]
        per-domain score cutoff:    [none]
        per-sequence Eval cutoff:   <= 10        
        per-domain Eval cutoff:     [none]
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        
        Query HMM:   smart_znf_c2h2_seed_alignment
        Accession:   [none]
        Description: [none]
          [No calibration for HMM; E-values are upper bounds]
        
        Scores for complete sequences (score includes all domains):
        Sequence                           Description          Score    E-value  N 
        --------                           -----------          -----    ------- ---
        chr14:0:125194864:frame=2:strand=0                     1299.4          0 301
        chr14:0:125194864:frame=0:strand=0                     1276.0          0 296
        chr14:0:125194864:frame=1:strand=0                     1093.6          0 306
        chr14:0:125194864:frame=1:strand=1                     1052.7          0 281
        chr14:0:125194864:frame=0:strand=1                      818.8     2e-246 239
        chr14:0:125194864:frame=2:strand=1                      775.9   1.7e-233 254
        
        Parsed for domains:
        Sequence                           Domain  seq-f seq-t    hmm-f hmm-t      score  E-value
        --------                           ------- ----- -----    ----- -----      -----  -------
        chr14:0:125194864:frame=1:strand=1 155/281 23493154 23493176 ..     1    23 []    25.8  1.1e-07
        
        Alignments of top-scoring domains:
        chr14:0:125194864:frame=1:strand=1: domain 155 of 281, from 23493154 to 23493176: score 25.8, E = 1.1e-07
                           *->ykCpigCgksfssk.aLkrHmrvH<-*
                              ++C i C +sfs+ ++L+ H+r+H   
          chr14:0:12 23493154    FQCRI-CMRSFSRSdHLTTHIRTH    23493176
        
        
        Histogram of all scores:
        score    obs    exp  (one = represents 1 sequences)
        -----    ---    ---
        > 775      6      -|======                                                     
        
        
        % No statistical fit available
        
        Total sequences searched: 6
        
        Whole sequence top hits:
        tophits_s report:
             Total hits:           6
             Satisfying E cutoff:  6
             Total memory:         21K
        
        Domain top hits:
        tophits_s report:
             Total hits:           1677
             Satisfying E cutoff:  1677
             Total memory:         641K

        """
        test_str = [s + "\n" for s in test_str.replace("\n        ", "\n").split("\n")]
        res = list(parse_hmm_search_results (test_str))
        print json.dumps(res, indent = 4)

    def test_function3(self):
        """
            test 
        """
        test_str = """hmmsearch - search a sequence database with a profile HMM
        HMMER 2.3.2 (Oct 2003)
        Copyright (C) 1992-2003 HHMI/Washington University School of Medicine
        Freely distributed under the GNU General Public License (GPL)
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        HMM file:                   /net/cpp-group/Leo/inprogress/6frame_translation_hmm/smart_znf_c2h2_seed_alignment.hmm [smart_znf_c2h2_seed_alignment]
        Sequence database:          /net/cpp-group/Leo/inprogress/6frame_translation_hmm/working_dir/mm_mm9_softmasked/mm_mm9_softmasked.00021.chr14.translated_fa
        per-sequence score cutoff:  [none]
        per-domain score cutoff:    [none]
        per-sequence Eval cutoff:   <= 10        
        per-domain Eval cutoff:     [none]
        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        
        Query HMM:   smart_znf_c2h2_seed_alignment
        Accession:   [none]
        Description: [none]
          [No calibration for HMM; E-values are upper bounds]
        
        Scores for complete sequences (score includes all domains):
        Sequence                           Description          Score    E-value  N 
        --------                           -----------          -----    ------- ---
        chr3:299990000:227962102:frame=2:strand=0 313/1082 8018360 8018382 ..     1    23 []    16.9  4.9e-05
        chr14:0:125194864:frame=2:strand=1                      775.9   1.7e-233 254
        
        Parsed for domains:
        Sequence                           Domain  seq-f seq-t    hmm-f hmm-t      score  E-value
        --------                           ------- ----- -----    ----- -----      -----  -------
        chr3:299990000:227962102:frame=2:strand=0 313/1082 8018360 8018382 ..     1    23 []    16.9  4.9e-05
        chr14:0:125194864:frame=1:strand=1 155/281 23493154 23493176 ..     1    23 []    25.8  1.1e-07
        
        Alignments of top-scoring domains:
        chr14:0:125194864:frame=1:strand=1: domain 155 of 281, from 23493154 to 23493176: score 25.8, E = 1.1e-07
                           *->ykCpigCgksfssk.aLkrHmrvH<-*
                              ++C i C +sfs+ ++L+ H+r+H   
          chr14:0:12 23493154    FQCRI-CMRSFSRSdHLTTHIRTH    23493176
        
        chr3:299990000:227962102:frame=2:strand=0: domain 313 of 1082, from 8018360 to 8018382: score 16.9, E = 4.9e-05
                           *->ykCpigCgksfssk.aLkrHmrvH<-*
                              y+ ++ Cgk f+ + +L +H r+H
          chr3:29999 8018360    YEYKH-CGKPFTERcSLFQHERIH    8018382
        
        Histogram of all scores:
        score    obs    exp  (one = represents 1 sequences)
        -----    ---    ---
        > 775      6      -|======                                                     
        
        
        % No statistical fit available
        
        Total sequences searched: 6
        
        Whole sequence top hits:
        tophits_s report:
             Total hits:           6
             Satisfying E cutoff:  6
             Total memory:         21K
        
        Domain top hits:
        tophits_s report:
             Total hits:           1677
             Satisfying E cutoff:  1677
             Total memory:         641K
        
        """
        test_str = [s + "\n" for s in test_str.replace("\n        ", "\n").split("\n")]
        res = list(parse_hmm_search_results (test_str))
        print json.dumps(res, indent = 4), "!!"
        res = self.adjust_results_with_contig_offsets(res)
        print json.dumps(res, indent = 4), "!!!"

#
#   debug code not run if called as a module
#     
if __name__ == '__main__':
    if sys.argv.count("--debug"):
        sys.argv.remove("--debug")
    unittest.main()



    
